version 14
clear all
adopath + ../../../lib/stata/gslab_misc/ado
adopath + ../../../lib/stata/ado
set more off
loadglob using  "variables_for_tables.txt" // FILE CONTAINS GLOBAL VARS. USED IN THIS CODE

program main_predict
    clean_data_predict, dep_var(enrollee_9m)
	gen one = 1
	classification, indep_var(age_80plus non_english                			///
                                   city_pitt male 								///
                                   race_white race_black race_others            ///
                                   Ibefore_2011                                 ///
                                   spending hospital SNF ED doctor              ///
                                   Ifull_year chronic one)
end

program clean_data_predict
    syntax, dep_var(str)
    use ../temp/exhibit_analysis, clear
	cap drop Itrain Ipredict Itemp xb xbd threshold 
    cap drop predict_enroll
	cap drop one 
	cap drop logit_* id_* prob_* rmax id count			
    cap drop *_xb
    gen Itrain = control == 1
    gen Ipredict = SP == 1
             
    save ../temp/predict_enroll, replace
end

program classification
	syntax, indep_var(str)
    foreach v in `indep_var' {
        assert mi(`var') == 0
    }

	qui gen Itemp = enrollee_9m == 1 & Itrain == 1
	qui replace Itemp = . if Itrain != 1
	qui count if Itrain == 1
	local N = r(N)
	qui count if Itemp == 1
	local s = log( (r(N)/`N'+ 1e-3) / (1 - r(N)/`N') )
	di "S is `s'"
		
	logit Itemp `indep_var' if Itrain == 1, noconstant
	predict xb, xb
 
	foreach var in `indep_var' {
		qui local `var'_b = _b[`var']
		di "Coef of `var' is ``var'_b'"
		qui gen `var'_xb = ``var'_b'*`var'
		matrix coef = nullmat(coef) \ round(``var'_b', 0.0001)
	}
	
	qui count if enrollee_9m == 1 & Itrain == 1
	local goal = r(N)
	qui count if xb > `s' & Itrain == 1
	di "in prediction `r(N)' obs Ipredict = 1"
	
	sum xb if Itrain == 1
	gen threshold = `s' - r(mean)
	gen xbd = xb + threshold
	
	gen predict_enroll = .
    qui replace predict_enroll = xbd > `s'
    
	forval it = 1/100 {
        di "it `it'"
		qui count if predict_enroll == 1 & Itrain == 1
		di "Now is `r(N)', goal is `goal'"
		qui replace threshold = threshold + 0.05 * (`s'  - log( (r(N)/`N'+ 1e-3) / (1 - r(N)/`N') ) )
		qui unique threshold
		assert r(sum) == 1
		qui replace xbd = xb + threshold

		qui cap drop temp
		qui replace predict_enroll = xbd > `s'
	}
	
    replace predict_enroll = . if Ipredict != 1
    save ../temp/predict_enroll, replace
end
	
main_predict
